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Abstract. We present spin-independent and spin-spin interquark potentials for charmonium states, that are calculated using a 
relativistic heavy quark action for charm quarks on the PACS-CS gauge configurations generated with the Iwasaki gauge action 
and 2-1-1 flavors of Wilson clover quark. The interquark potential with finite quark masses is defined through the equal-time 
Bethe-Salpeter amplitude. The light and strange quark masses are close to the physical point where the pion mass corresponds 
to ~ 156(7) MeV, and charm quark mass is tuned to reproduce the experimental values of r\c and J/\^f states. Our 
simulations are performed with a lattice cutoff of ~ 2.2 GeV and a spatial volume of (3 fm)^. We solve the nonrelativistic 
Schrodinger equation with resulting charmonium potentials as theoretical inputs. The resultant charmonium spectrum below 
the open charm threshold shows a fairly good agreement with experimental data of well-established charmonium states. 
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INTRODUCTION 

The heavy-quark (2)-antiquark {Q) potential is an important quantity to understand many properties of the heavy 
quarkonium states. The dynamics of heavy quarks can be described well within the framework of nonrelativistic 
quantum mechanics, because of their masses being much larger than the QCD scale (Aqcd)- Indeed the constituent 
quark potential models with a QCD-motivated QQ potential have successfully reproduced the heavy quarkonium 
spectra and also decay rates below open thresholds [1, 2, 3]. 

In the nonrelativistic potential (NRp) models, the heavy quarkonium states such as charmonium and bottomnium 
are well understood as is a quark-antiquark pair bound by the Coulombic induced by perturbative one-gluon exchange, 
plus linearly rising potential. The former dominates in short range, while the latter describes the phenomenology of 
confining quark interactions at large distances [1]. This potential is called as Cornell potential and its functional form 
is given by V (r) = — j ^ + ar -f Vb where and a denote the strong coupling constant and the string tension, and Vb 
is the constant term associated with a self-energy contribution of the color sources. In the NRp models, spin-dependent 
potentials are induced as relativistic corrections in powers of the relative velocity of quarks, and their functional forms 
are also determined on the basis of perturbative one-gluon exchange as the Fermi-Breit type potential [4]. However 
the validity of the phenomenological spin-dependent potentials determined within the perturbative method would be 
limited only at short distances and also in the vicinity of the heavy quark mass limit. This may cause large uncertainties 
in the predictions for higher heavy quarkonium states obtained in the NRp models. 

Lattice QCD simulations offer a strong tool to understand the properties of QQ interactions. Indeed, both the static 
QQ potential and its corrections of order ^(1 /m^) as the spin dependent potentials have been precisely determined 
from Wilson loops using lattice QCD simulations with the multilevel algorithm [5, 6]. Although the lattice QCD 
calculations within the Wilson loop formalism support a shape of the Cornell potential [7], the leading spin-spin 
potential determined at ^(l/m^) gives an attractive interaction for the higher spin states [8, 9], in contradiction with 
a repulsive one that is demanded by phenomenological analysis. The higher order corrections beyond the next-to- 
leading order are required to correctly describe the conventional charmonium spectrum, because the inverse of the 
charm quark mass would be far outside the validity region of the 1 /mg expansion [10]. In addition, practically, the 
multilevel algorithm is quite difficult to be implemented in dynamical lattice QCD simulations. 

Under this situation, we employ the new method proposed in our previous works [10, 11, 12] in order to obtain 
the proper interquark potentials fitted in the NRp models. The interquark potential and also the quark kinetic mass 
are defined by the equal-time and Coulomb gauge Bethe-Salpeter (BS) amplitude through an effective Schrodinger 
equation. This new method enables us to determine the interquark potentials including spin-dependent terms al finite 
quark masses from first principles of QCD, and then can fix all parameters that are needed in the NRp models. 
Furthermore, there is no restriction to extend to the dynamical calculations. Hereafter we call the new method as BS 



amplitude method. 

Once we obtain the reliable QQ potentials from lattice QCD, we can solve the nonrelativistic Schrodinger equa¬ 
tion with “lattice-determined potentials” as theoretical inputs, and obtain many physical observables such as mass 
spectrum. In this proceedings, we present not only the charmonium potentials calculated with almost physical quark 
masses using the 2-1-1 flavor PACS-CS gauge configurations [13], but also the resultant charmonium mass spectrum 
computed from the NRp model with the lattice-determined potentials, where there are no free parameters including the 
Vo and quark mass. The simulated pion mass « 156(7) MeV is almost physical. For the heavy quarks, we employ 
the relativistic heavy quark action that can control large discretization errors introduced by large quark mass [14]. 


FORMALISM 


In this section, we briefly review the BS amplitude method utilized to calculate the interquark potential with the finite 
quark mass. This is an application based on the approach originally used for studing the hadron-hadron potential, which 
is defined through the equal-time BS amplitude [15, 16]. More details of determination of the interquark potential are 
given in Ref. [10]. 

In lattice simulations, we measure the following equal-time QQ BS amplitude in the Coulomb gauge for the 
quarkonium states [17, 18]: 

.()r(r) =^(0|e(x)re(x + r)|ee;7^C), (1) 

X 

where r is the relative coordinate between quark and antiquark at time slice t. The Dirac ymatrices F in Eq. (1) specify 
the spin and the parity of meson states. For instance, 75 and 7 correspond to the pseudoscalar (PS) and the vector (V) 
channels with = 0 ^ and = 1 , respectively. A summation over spatial coordinates x projects onto zero 

total momentum. The r-dependent amplitude, ^(r), is here called BS wave function. The BS wave function can be 
extracted from four-point correlation function at large time separation. Also, the corresponding meson masses Mr can 
be read off from the asymptotic large-time behavior of two-point correlation functions. In this proceedings, we focus 
only on the 5-wave charmonium states ( 77 ^ and obtained by appropriate projection to the representation in 

cubic group [19]. 

The BS wave function satisfies an effective Schrodinger equation with a nonlocal and energy-independent interquark 
potential U [15, 20, 21] 

J dr'U(ry)(j)rir')=Er(j)r{r), ( 2 ) 

where 71 is the reduced mass of the QQ system. The energy eigenvalue of the stationary Schrodinger equation is 
supposed to be Mp — 2mg. If the relative quark velocity v = | V/mg] is small as v ^ 1, the nonlocal potential U can 
generally expand in terms of the velocity v as t/(r',r) = {V(r) -|-Vs(r)Sg-Sg-|-VT(r) 5 i 2 -l-VLs(r)L-S-|- ^(v^)}5(r' — 
r) where 5 i 2 = (Sg • f)(Sg • f) — Sg • Sg/3 with f = r/r, S = Sg-I-Sg and L = r x (—/V) [15]. Here, V, Vs, Vt and Vps 
represent the spin-independent central, spin-spin, tensor and spin-orbit potentials, respectively. 

The Schrodinger equation for 5-wave is simplified as 

{ - ^ + V (r) + Sg • SgVs (r)} ())r(r) = £r0r(r) (3) 

at the leading order of the v-expansion. Here, we essentially follow the NRp models, where the J/xj/ state is purely 
composed of the 15 wave function. 

The spin operator Sg • Sg can be easily replaced by expectation values —3/4 and 1 /4 for the PS and V channels, 
respectively. Then, the spin-independent and spin-spin QQ potentials can be evaluated through the following linear 
combinations of Eq.(3): 
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(4) 
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where Eave = 47ave ~ 2mg and Ehyp = My — Mps. The mass Mave denotes the spin-averaged mass as ^Mps -f ^My. The 
derivative is defined by the discrete Laplacian. 



TABLE 1. Parameters of 2 +1 -flavor dynamical QCD gauge field configurations generated by PACS-CS collab¬ 
oration [13]. The columns list number of flavors, lattice volume, the j3 value, hopping parameters (light, strange), 
approximate lattice spacing (lattice cut-off), spatial physical volume, pion mass, number of configurations to be 
analyzed. 


Nf 

L^xT 

P 

Kud 

Ks 

a [fm] (a * 

[GeV]) 

La [fm] 

Mn [MeV] 

# configs. 

2+1 

32^ X 64 

1.9 

0.13781 

0.13640 

0.0907 Q 

e 2.176) 

^2.90 

s:il56 

198 


TABLE 2. The hopping parameter Kq and RHQ 
parameters used for the charm quark. 


Kc 

V 

rs 

Cb 

Ce 

0.10819 

1.2153 

1.2131 

2.0268 

1.7911 


The kinetic quark mass is an important quantity in the determination of the interquark potentials since Eqs. (4) 
and (5) require an information of the kinetic quark mass mg. In our previous work [10, 11, 12], we propose to 
calculate the quark kinetic mass through the large-distance behavior in the spin-spin potential with the help of the 
measured hyperfine splitting energy of 15 states in heavy quarkonia. Under a simple, but reasonable assumption as 
limr^ooV 5 (r) = 0 which implies there is no long-range correlation and no irrelevant constant term in the spin-spin 
potential, Eq. (5) is rewritten as 


_ -1 V^^s(r)\ 

r^Ehypl 0v(O 0ps(?-) j 


( 6 ) 


and then we can estimate the kinetic quark mass from asymptotic behavior of Eq. (6) in long range region. 


LATTICE SETUP 

The computation of the charmonium potential in this study is performed on a lattice T = 31? X 64 using the 2+1 
flavor PACS-CS gauge configurations [13] generated by non-perturbatively ^(a)-improved Wilson quark action with 
csw = 1.715 [22] and Iwasaki gauge action at P — 1.90 [23], which corresponds to a lattice cutoff of = 2.176(31) 
GeV (a = 0.0907(13)fm). The spatial lattice size then corresponds to La « 3 fm. The hopping parameters for the light 
sea quarks ( K'„^,k;s}={ 0.13781, 0.13640} provide = 156(7) MeV andM/f = 554(2) MeV [13]. Table 1 summarizes 
simulation parameters of dynamical QCD simulations used in this work. Although the light sea quark masse is slightly 
off the physical point, the systematic uncertainty due to this fact could be extremely small in this project. Our results 
are analyzed on all 198 gauge configurations. All gauge configurations are fixed to Coulomb gauge. 

In order to control discretization errors induced by large quark mass, we employ the relativistic heavy quark (RHQ) 
action [14] that removes main errors of ^{\p\a), ^((mofl)”) and ^?{\p\a{moa)") from on-shell Green’s functions. The 
RHQ action is the anisotropic version of the i^{a) improved Wilson action with five parameters v, cb and ce, 
called RHQ parameters (for more details see Ref. [14, 24]). The RHQ action utilized here is a variant of the Eermilab 
approach [25] (See also Ref. [26]). 

The parameters cb and ce in RHQ action are determined by tadpole improved one-loop perturbation theory [24]. 
Eor V, we use a nonperturbatively determined value, which is adjusted by reproducing the effective speed of light Ceff 
to be unity in the dispersion relation ^^(p^) = +Cgg|pp for the spin-averaged 15-charmonium state, since the 
parameter v is sensitive to the size of hyperfine splitting energy [27]. We choose jq to reproduce the experimental 
spin-averaged mass of 15-charmonium states Ml^g{lS) = 3.0678(3) GeV. To calibrate adequate RHQ parameters, 
we employ a gauge invariant Gauss smearing source for the standard two-point correlation function with four finite 
momenta. As a result, the relevant speed of light in the dispersion relation is consistent with unity within statistical 
error: = 1.04(5). Our chosen RHQ parameters are summarized in Table 2. 

Using tuned RHQ parameters, we compute two valence quark propagators with wall sources located at different time 
slices fs/a = 6 and 57 to increase statistics. Two sets of two and four-point correlation functions are constructed from 
the corresponding quark propagators, and folded together to create the single correlation function. Dirichlet boundary 
condition is imposed for the time direction to eliminate unwanted contributions across time boundaries. 



TABLE 3. Masses of low-lying charmonium states 
calculated from two-point functions, the spin-averaged 
mass and hyperfine splitting energy of IS charmonium 
states. The fitting ranges and values of ;f^/d.o.f. are 
also included. Results are shown in units of GeV. 


state (7^*') 

ht range 

mass [GeV] 

;t:^/d.o.f. 

ric (0-+) 

[33:47] 

2.9851(5) 

0.70 

7/vr(l-+) 

[33:47] 

3.0985(11) 

0.62 

XcO (0++) 

[14:26] 

3.3928(59) 

0.66 

Xcl (1++) 

[14:26] 

3.4845(62) 

1.03 

he(l + -) 

[14:26] 

3.5059(62) 

0.63 

Mave(15) 

- 

3.0701(9) 

- 

£hyp(15) 

- 

0.1138(8) 

- 



FIGURE 1. The reduced QQ BS wave functions of the rjc (circles) and J/lj/ (squares) states, shown as a function of the spatial 
distance r. The data points are taken along r vectors which are multiples of three directions (1,0,0), (1,1,0) and (1,1,1). 


Low-lying charmonium masses of rjc, J/yf, he, XcO and Xci are obtained by weighted average of the effective mass 
in the appropriate range. The effective mass is dehned as 


Mr{t)= log 


Gr{t,ts) 
Gr{t + l,fs) ’ 


(7) 


where Gr(f,fs) is a two-point function obtained by setting r to be zero in four-point function Gr(r,f,fs). In Table 2, 
we summarize resultant charmonium masses together with ht ranges used in the hts and j^/d.o.f. values. We take 
into account a correlation between effective masses measured at various time slices in the ht. The statistical errors are 
estimated by the jackknife method. 

Low-lying charmonium masses calculated in this study below DD threshold are all close to the experimental 
values, though the hyperhne mass splitting Mhyp = 0.1124(9) GeV is slightly smaller than the experimental value, 
= 0.1166(12) GeV [28]. Note that here we simply neglect the disconnected diagrams in two-point correlation 
functions. The several numerical studies reported that the contributions of charm annihilation to the hyperhne splitting 
of the 15-charmonium state are sufficiently small, as of order 1 —4 MeV. [29, 30, 31], 


DETERMINATION OF INTERQUARK POTENTIAL 
QQ BS wave function 

Fig. 1 shows the QQ BS wave functions of IS charmonium states (rje and J/y/ states). The BS wave functions 
are dehned by Eq.(l) and normalized as ~ 1- reduced wave function ur{r) for displaying the wave 

function: Mr(^) = Practically we take average of the BS wave function by weight over time slices 33 < f/a < 47 








FIGURE 2. The determination of quark kinetic mass within the BS amplitude method. The values of — (V^(|v/i|V ~ 
V^#/#)/Ehyp as a function of the spatial distance r are shown in this figure. The quark kinetic mass niQ is obtained from 
the long-distance asymptotic values of — (V^(/v/(^ — V^0p/0p)/£’l,yp. Horizontal solid line indicates a value of quark kinetic mass 
obtained by fitting a asymptotic constant in the range 0.54 fm < r < 1.10 fm. A shaded band indicates a statistical error estimated 
by estimated by jackknife method. 


where effective mass plots for 15-charmonium states show plateaus and excited state contaminations are expected to 
be negligible. In Fig. 1 we display data points of Mr(r) calculated at r vectors which are multiples of (1,0,0), (1,1,0) 
and (1,1,1). Hereafter we focus on lattice data taken in three directions for any quantities. 

We find that a sign of rotational symmetry breading found in the QQ BS wave functions is sufficiently small in our 
calculation. The resulting wave functions become isotropic with the help of a projection to the sector of the cubic 
group that corresponds to the 5-wave in the continuum theory (Fig. 1). 


quark kinetic mass 


In our formalism, the kinetic mass of the charm quark is determined self-consistently within the BS amplitude 
method as well [11]. The quark kinetic mass defined in Eq. (6) is calculated from asymptotic behavior of the quantity 
— {V^ (j)Y/(j)y — (j)p/(j)p) / E^yp at long distances. Fig. 2 illustrates the determination of quark kinetic mass iuq for the 
charmonium system. 

For the derivative, we use the discrete Laplacian operator defined in polar coordinates as 




2 ^{r + d)- ^(r-d) ^ ^{r + d) + ^{r-d)-2^ 
r 2d d^ 


where r is the absolute value of the relative distance as r = |r | and a is a spacing between grid points along differentiate 
directions. In the on-axis (r (1,0,0)) and the two off-axis directions (r (1,1,0) and (1,1,1)), the effective grid 
spacings correspond to d = a, v/2a, y/Sa, respectively. 

The differences of ratios V^0r/())r at each r are obtained by a constant fit to the lattice data with a reasonable 
X^/d.o.f. value over the range of time slices where two-point functions exhibit the plateau behavior (33 <tla< 47). 
Then the values of mq are determined for each directions from asymptotic values of —(V^tjtvl(j>v — V^<)>/>/<)>p)/Ehyp in 
the range of 6 < r/a < 7v/3 where V’s(r) should vanish. Finally we average them over three directions, and then obtain 
iriQ — 1.784(23) (6) (20) GeV. The first error is statistical, given by the jackknife analysis. In the second error, we quote 
a systematic uncertainty due to rotational symmetry breaking by taking the largest difference between the average value 
and individual ones obtained for specific directions. The third one represents the systematic uncertainties due to choice 
of fmin of the time range used in the fits. We vary fn,in over range 33 — 41 and then quote the largest difference from the 
preferred determination of mg. 
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FIGURE 3. Central spin-independent and spin-spin charmonium potentials calculated from the BS wave functions in the 
dynamical QCD simulation with almost physical quark masses. In the upper panel, we show the spin-independent potential V(r). A 
solid (dot-dashed) curve is the fit results with the Cornell (Cornell plus log) form. The shaded bands show statistical uncertainties in 
the fitting procedure where the jackknife analysis is used. Note that the spin-averaged eigen-energy of 15-charmonium state £ave is 
not subtracted in this figure. A horizontal line indicates the level of open-charm {D^D^) threshold ~ 3729 MeV. In the lower panel, 
we show the spin-spin potential Vs{r). A solid (dot-dashed) curve corresponds to fitting results with exponential (Yukawa) form. 
The inset shows a magnified view. In both plots, the phenomenological potentials adopted in a NRp model [3] are also included as 
dashed curves for comparison. 


TABLE 4. Summary of the Cornell parameters and the 
quark mass determined by the BS amplitude method. For 
comparison, ones adopted in a phenomenological NRp 
model [3] and ones of the static potential obtained from 
Polyakov line correlations are also included. In the first col¬ 
umn, the quoted errors indicate the sum of the statistical and 
systematic added in quadrature. 



This work 

Polyakov lines 

NRp model 

A 

0.713(83) 

0.476(81) 

0.7281 

x/ct [GeV] 

0.402(15) 

0.448(16) 

0.3775 

mg [GeV] 

1.784(31) 

OO 

1.4794 


Spin-independent interquark potential 

Once the quark kinetic mass is determined, we can easily calculate the central spin-independent and spin-spin 
charmonium potentials from the QQ BS wave function through Eqs. (4) and (5). First, we show a result of the spin- 
independent charmonium potential V (r) in Fig. 3. The constant energy shift Fave is not subtracted. At each distance r, 
the values of interquark potentials V (r) and Vs{r) are practically determined by constant fits to data points over time 



















slices where two-point functions exhibit the plateau behavior. The correlations between data points at different time 
slices are taken into account in the fitting process. 

The charmonium potential calculated by the BS amplitude method from dynamical attice QCD simulations properly 
exhibits the linearly rising potential at large distances and the Coulomb-like potential at short distances. The finite mg 
corrections could be encoded into the Cornell parameters, although the charm quark mass region would be beyond 
the radius of convergence for the systematic 1 /mg expansion. Therefore, as first step, we simply adopt the Cornell 
parametrization to fit the data of the spin-independent central potential: V{r) = —j: + (Jr + Vo with the Coulombic 
coefficient A, the string tension a, and a constant Vq. 

All fits are performed individually for each three directions over the range [rniin/a,?‘max/a] = [4 : 7-\/3]. We 
minimize the j^/d.o.f including the covariance matrix. Resulting Cornell parameters of the charmonium potential 
are A = 0.713(26)(38)(31)(62) and ^/a — 0.402(6)(4)(9)(9) MeV with j^/d.o.f. « 3.2. The first error is statistical 
and the second, third and forth ones are systematic uncertainties due to the choice of the differentiate direction, 
fmin ^min^ respectively. The resulting Cornell parameters are summarized in Table 4. Also we include both 
phenomenological ones adopted in the NRp model [3] and of the static potential obtained from Polyakov loop 
correlations. The latter is calculated using the same method as in Ref. [13]. Additionally we calculate the Sommer 
parameter defined as ro = ■y/(1.65 —A)/o, and then obtain ro = 0.476(6)(11)(3)(6) fm, which is fairly consistent 
with the value quoted in Ref. [13]. 

As shown in Table 4, a gap for the Cornell parameters between the conventional static potential from Wilson- 
loops (Polyakov-loops) and the phenomenological potential used in the NRp models seems to be filled by our new 
approach, which nonperturbatively accounts for a finite quark mass effect. In the charmonium potential from the BS 
wave function, a Coulomb-like behavior is enhanced and the linearly rising force is slightly reduced due to finite charm 
quark mass effects. For the spin-independent central interquark potential, the 1 /mg expansion within the Wilson-loop 
approach converges in the heavy quark mass region of mg >1.8 GeV. Indeed, as reported in Ref. [32], the static QQ 
potential and its 1 jriiQ corrections calculated in Ref. [33] agree with the charmonium potential obtained from the BS 
amplitude method. 

In order to provide a more adequate fit to the lattice data, we try to employ an alternative functional form adding a 
log term to the Cornell potential: 

y(r) =---|-(7r-f Vb-l-Blog(rA) (8) 

r 

where A is simply set to be lattice cutoff a^*. Such log term as 1 /mg corrections to the spin-independent potential is 
reported in Ref. [34]. Resulting parameters are A = 0.194(137) (33) (36) ( 66 ), x/a = 0.300(38) (19) (20) (21) GeV and 
B = 0.390(113)(20)(39)(61) GeV. with j^/d.o.f. sa 2.3. Fitting range is determined to minimized a j^/d.o.f. value 
taking into account the correlation, and then we choose rmax/a] = [3 : 7-\/3]. 

The finite quark mass corrections to spin-independent potential give only a minor modification in the NRp models. 
In the upper panel of Fig. 3 the solid (dot-dashed) curve is given by the fitting the data to the Cornell form (Cornell 
plus log form). The phenomenological potential used in the NRp models [3] is also plotted as a dashed curve for 
comparison. The charmonium potential obtained from lattice QCD is similar to the one used in the NRp models, 
although a slope of the charmonium potential in the long range is barely larger than the phenomenological one. 

It is worth mentioning that a string breaking-like behavior found in the range r < 1.1 fm is unreliable. In principle, 
string breaking due to the presence of dynamical quarks is likely to be observed. The signal-to-noise ratio however 
becomes worse rapidly for the spin-independent potential as spatial distance r increase because of the localized wave 
function. The lattice data of the potential near the spatial boundary are also sensitive to finite volume effects. Therefore, 
at least, calculations of the higher charmonium near the open charm threshold using a larger lattice are required for 
observing the string breaking. Their wave functions are extended until the string breaking sets in. 


Spin-Spin potential 

The lower plot of Fig. 3 shows the spin-spin charmonium potential obtained from the BS amplitude method with 
almost physical quark masses. The spin-spin potential exhibits the short-range repulsive interaction, which is required 
to leads heavier mass to the higher spin state in hyperfine multiplets. In contrast of the case of the spin-independent 
potential, the spin-spin potential obtained from BS wavefunction is absolutely different from a repulsive 5-function 
potential generated by perturbative one-gluon exchange [4]. Such contact form oc 5 (r) of the Fermi-Breit type potential 
is widely adopted in the NRp models [2]. 



TABLE 5. Results of fitted parameters for the spin-spin poten¬ 
tial with the exponential and Yukawa forms. The quoted errors 
are statistical only. In the case of the spin-spin potential, we use 
only on-axis data. 


Functional form 

a 

j8 

X^/d.o.i. 

Exponential 

2.15(7) GeV 

2.93(3) GeV 

2.0 

Yukawa 

0.815(27) 

1.97(3) GeV 

1.7 


The QQ interaction is not entirely due to one-gluon exchange so that spin-spin potential is not necessary to be a 
simple contact form 0 = 5(r). Indeed, the finite-range spin-spin potential described by the Gaussian form is adopted 
by the phenomenological NRp model in Ref. [3], where many properties of conventional charmonium states at 
higher masses are predicted. This phenomenological spin-spin potential is also plotted in the lower plot of Fig. 3 
for comparison. There is a slight difference at very short distances, although the range of spin-spin potential calculated 
from the BS amplitude method is similar to the phenomenological one. 

To examine an appropriate functional form for the spin-spin potential, we try to fit the data with several functional 
forms, and explore which functional form can give a reasonable fit over the range of rja from 2 to 7 a/3. As a results, 
the long-range screening observed in the spin-spin potential is accommodated by the exponential form or the Yukawa 
form: 



aexp(—j3r) 
aexp(—j3r)/r 


Exponential form 
Yukawa form 


(9) 


All results of correlated fits are summarized in Table 5. We also try to fit the data with the Gaussian form that is 
often employed in the NRp models, however it provides an unreasonable j^/d.o.f. value. Note that we here use only 
the on-axis data which are expected to less suffer from both the rotational symmetry breaking and discretization error, 
because fit results obtained in each direction significantly disagree with each other. We need the finer lattice to make a 
solid conclusion regarding the shape of the spin-spin potential and also systematic uncertainties due to the rotational 
symmetry breaking. 


NONRELATIVISTIC POTENTIAL MODEL WITH LATTICE INPUTS 

Using the quark kinetic mass and the charmonium potentials determined by first principles of QCD, we can solve 
the nonrelativistic Schrodinger equation for the bound cc systems as same as calculations in the NRp models. In the 
BS amplitude method, a value of the difference Vo — Eave is directly obtained as the constant term in spin-independent 
charmonium potential, while the value of Eave is calculated through £ave = — ^nig. However statistical uncertainty 

of niQ is somewhat large compared to an error of Vb —Eave: here these are £ave = 0.508(69) GeV, mg = 1.789(34) GeV 
and Vo — £ave = —0.146(13) GeV. To reduce statistical uncertainties, we therefore solve the following Schrodinger 
equation shifted by a constant energy —£ave: 

I ~ 2 ^ + ^sui^) \ '^suir) = E'^ijUsuir) (10) 

( mg mgr^ J 

where Vgj^{r) =VsLj{f) — ^ave andE^^^ = £ave- The interquark potential depends on the channel of charmonium 

states with S, L and J. Desired charmonium masses are obtained by merely adding to the spin-averaged mass Mave 
which is obtained from the standard lattice spectroscopy with high accuracy: M^u = Mave +E 5 /J = '2-mq + E^u. 

The resulting potentials from lattice QCD are discretized in space [35]. Therefore, instead of solving continuum-type 
Schrodinger equation, we practically solve eigenvalue problems as 


V/ 2-1 




= Eu„ 


n=\ 


( 11 ) 
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FIGURE 4. Mass spectrum of charmonium states below and near the open-charm threshold. The vertical scale is in units of 
MeV. Labels of [L\j ) are displayed in lower (upper) horizontal axis. Rectangular boxes indicate predictions from the NRp 
models with theoretical inputs based on lattice QCD and their errors which are the sum of the statistical and systematic added in 
quadrature. Solid lines indicate experimental values of well established charmonium states, while square symbols represent results 
of the standard lattice spectroscopy. A horizontal line shows the open-charm threshold. A symbol of ^Pj denotes the spin-weighted 
average of spin-triplet ^Pj states as Mjj = )/9- 


with a symmetric matrix defined in one of three specific directions 


Hnjt — 


1 


cf-mq 


2 + 


L[L+\) 


+ V'{no) 


Hn±l,n — 


( 12 ) 

(13) 


The boundary condition to the reduced wave functions = u(na) is simply set to mq = 0 and MyVs /2 = 0- In this work, we 
separately solve Eq. (11) in the directions of vectors r which are multiples of (1,0,0), (1,1,0) and (1,1,1). We prefer to 
use mainly on-axis data which is expected to receive smallest discretization etTors and systematic uncertainties due to 
rotational symmetry breaking, and quote the largest difference between on-axis and off-axis results as the systematic 
etTor due to the choice of direction. While statistical errors are estimated by the jackknife method. A systematic 
uncertainty stemming from the choice of time window is relatively small compared with other errors. Alternatively 
we can solve the Schrodinger equation in continuum space with the parameterized charmonium potential by empirical 
functional forms. This procedure however highly depends on choice of functional forms especially at short distances, 
and give additional systematic uncertainties to resultant spectrum. 

Fig. 4 shows the mass spectrum of the charmonia below 4200 MeV. Theoretical spectrum plotted as rectangular 
boxes are given by solving the discrete nonrelativistic Schrodinger equations with theoretical inputs. Quoted errors of 
charmonium masses are statistical and systematic uncertainties combined in quadrature. For purpose of comparison, 
both experimental values and results of the standard lattice spectroscopy are plotted together. The experimental values 
are taken from Particle Data Group [28]. At first glance, we find that theoretical calculations from the NRp model 
with lattice inputs are in fairly good agreement with not only the lattice spectroscopy, but also experiments below 






















TABLE 6. Charmonium mass spectrum is summarized in units of MeV. 
The labels of AVE and HYP in a column of “state” for S-states denotes 
the spin-averaged mass (Mi + SMs) /4 and hyperfine splitting mass 
Ml— Ms . Experimental data (denoted as Exp.) are taken from Particle 
Data Group, rounded to 1 MeV [28]. There are two lattice QCD results. 
Eirst one is given by the usual spectroscopy, and second one is a result 
calculated by solving the Schrbdinger equation with the charmonium 
potentials determined from lattice QCD. For the second, first error is 
statistical, and second error is systematic error due to rotational symmetry 
breaking. For spin triplet states ^[L\j, the spin-weighted average (^[T]y) 
are also included. 


state 

Exp. 

Lattice QCD 

NRp model 



spectroscopy 

BS amplitude 



2981 

2985(1) 

2985(2)(1) 

2982 

y/v/(i3si) 

3097 

3099(1) 

3099(2)(1) 

3090 

AVE 

3068 

3070(9) 

3070(2)(1) 

3063 

HYP 

116 

114(1) 

113(1)(0) 

108 

r]c (2'So) 

3639 


3612(9)(7) 

3630 

V/(23Si) 

3686 


3653(12)(5) 

3672 

AVE 

3674 


3643(11)(5) 

3662 

HYP 

47 


41(6)(3) 

42 

r]c (3'So) 



4074(20)(70) 

4043 


4039 


4099(24)(98) 

4072 

AVE 



4092(22)(91) 

4065 

HYP 



25(15)(28) 

29 

he (l^Pi) 

3525 

3506(6) 

3496(7)(19) 

3516 

1(7 (l^Q) 

3525 


3503(7)(10) 

3524 

XcQ 

3415 

3393(6) 


3424 

Xc\ (l^A) 

3511 

3485(6) 


3505 

;fc2 

3556 



3556 

he {2^Pi) 



3927(16)(34) 

3934 

xn i2^Pj) 



3916(19)(31) 

3943 

XcO {2 ^Pq) 

3918 



3852 

Xel {2^Pi) 




3925 

Xe2 (2^Pl) 

3927 



3972 

riel (1‘D2) 



3783(12)(4) 

3799 

V7(l3Dy) 



3774(13)(2) 

3800 

Wil^Di) 

3773 



3785 

IX {1^02) 




3800 

Vr(l3D3) 




3806 

m (2'd2) 



4221(21)(72) 

4158 

W {2^Dj) 



4193(25)(88) 

4159 

¥(2^P>i) 

4153 



4142 

V {2^Dl) 




4158 

V {2^Di) 




4167 


open charm threshold. All results are also summarized in Table 6. In this study, we succeed in extracting only the 
spin-spin potential among spin-dependent interquark potentials. Thus at this stage we cannot predict the spin-orbit 
splitting which is led by the tensor and spin-orbit forces. In other wards, we can compute only the spin-averaged mass 
for excited states with higher angular momentum such as Xcj state. 

Our theoretical calculations for charmonium states below the open-charm threshold are in fairly good agreement 
with the experimental measurements. The point we wish to emphasize here is that our novel approach has no free 
parameters in solving the Schrodinger equation opposed to the phenomenological NRp model. All of the parameters 
appeared in the NRp model calculation are solely determined by lattice QCD simulations, where three light hadron 



masses (M-ji, Mk and Mq) are used for inputs to fix the lattice spacing and light quark hopping parameters. Only 
experimental values of rjc and J/xj/ masses in the charm sector are used to determine the charm quark parameters in 
the RHQ action. In this sense the new approach proposed here is distinctly different from the existing calculations with 
the phenomenological quark potential models. 


SUMMARY 

We have calculated the interquark potentials between charm quark and anti-charm quark almost on the physical 
point. The interquark potential at finite quark mass is defined through the equal-time Bethe-Salpeter wave function. 
Our simulations have been performed in the vicinity of the physical light quark masses, which corresponds to 
Alji — 156 MeV, using the PACS-CS gauge configurations generated with the Iwasaki gauge action and 2+\ flavors of 
Wilson clover quark. We use the relativistic charm quark tuned to reproduce the experimental values of T/y/ and Tj^ 
masses. The resulting spin-independent potential shows behavior of Coulomb plus linear form, and their parameters 
are close to values used in the traditional quark potential models. Also the string breaking due to existence of sea 
quarks is not observed. On the other hand, the spin-spin potential obtained from the dynamical simulations exhibits 
the short-range repulsive interaction. Its shape is quite different from the a repulsive 5-function potential induced by 
the one-gluon exchange which are usually adopted in the quark potential model. 

We have calculated the charmonium spectrum by solving nonrelativistic Schrodinger equation with the theoretical 
input of the spin-independent and spin-spin potentials and the quark kinetic mass. We simply solved the Schrodinger 
equation with Dirichlet boundary condition in a matrix manner. This approach enable us to directly use the raw data 
of the charmonium potential without introducing a phenomenological parameterization for the discretized potential 
data. We found an excellent agreement of low-lying charmonium masses between our results and the experimental 
data. We emphasize that our novel approach has no free parameters in solving the Schrodinger type equation opposed 
to conventional phenomenological quark potential models. As for inputs of lattice QCD, we essentially use three 
light hadron masses and Mq) for fixing the lattice spacing and light quark hopping parameters, and two 

charmonium masses (Mri^ and Mjj^) for determining the parameters of RHQ action. 

In order to precisely predict the mass spectrum above the open charm threshold, we must take into account the effects 
of not only the mass shift caused by mixing the QQ states with DD continuum, but also S-D mixing due to existence 
of the tensor force. However, in this work, we simply ignore these effects and also apply our new approach to the 
charmonium states above the open-charm threshold. The theoretical prediction of the nonrelativistic potential model 
with lattice inputs is basically consistent with the existing experimental data, although the systematic uncertainties due 
to the rotational symmetry breaking are rather large. For more comprehensive prediction including spin-orbit splittings, 
however, we must calculate all spin-dependent terms (spin-spin, tensor and spin-orbit forces). Especially the tensor 
force introducing the S-D mixing would shift even the masses of 15-states. Also the larger spatial extent is required to 
address the systematic uncertainties due to the finite size effect for the higher excited state that are supposed to possess 
wider wavefunction. 
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